nature astronomy 


Article https://doi.org/10.1038/s41550-023-02105-7 


Filamentary structures as the origin of blazar 
jet radio variability 


Antonio Fuentes ©'“, José L. Gomez ©‘), José M. Marti2*, Manel Perucho2*, 
Guang-Yao Zhao’, Rocco Lico®™, Andrei P. Lobanov ©”, Gabriele Bruni ©°, 
Yuri Y. Kovalev © ™, Andrew Chael ©, Kazunori Akiyama © °°", 

Katherine L. Bouman ©”, He Sun”, Ilje Cho ©', Efthalia Traianou', 

Teresa Toscano ©', Rohan Dahale ©’, Marianna Foschi®', 

Leonid I. Gurvits®“”, Svetlana Jorstad © ©”, Jae-Young Kim © *"2"°, 

Alan P. Marscher", Yosuke Mizuno ©7°"2, Eduardo Ros @° & 

Tuomas Savolainen”??? 


Received: 8 February 2023 


Accepted: 18 September 2023 


Published online: 26 October 2023 


® Check for updates 


Supermassive black holes at the centre of active galactic nuclei power some 
of the most luminous objects in the Universe. Typically, very-long-baseline 
interferometric observations of blazars have revealed only funnel-like 
morphologies with little information on the internal structure of the ejected 
plasma or have lacked the dynamic range to reconstruct the extended jet 
emission. Here we present microarcsecond-scale angular resolution images 
of the blazar 3C 279 obtained at 22 GHz with the space very-long-baseline 
interferometry mission RadioAstron, which allowed us to resolve the jet 
transversely and reveal several filaments produced by plasma instabilities 
in a kinetically dominated flow. The polarimetric properties derived 

from our high-angular-resolution and broad-dynamic-range images are 
consistent with the presence of a helical magnetic field threaded to the jet. 
We infer a clockwise rotation as seen in the direction of flow motion with an 
intrinsic helix pitch angle of -45° and a Lorentz factor of ~13 at the time of 
observation. We also propose a model to explain blazar jet radio variability 
in which emission features travelling down the jet may manifest as a result 
of differential Doppler boosting within the filaments, as opposed to the 
standard shock-in-jet model. Characterizing such variability is particularly 
important given the relevance of blazar physics from cosmic particle 
acceleration to standard candles in cosmology. 


We observed 3C 279 on 10 March 2014 at 22 GHz (1.3 cm) with the 
space very-long-baseline interferometry (VLBI) mission RadioAstron', 
a10 m space radio telescope onboard the Spektr-R satellite and an 
array of 23 ground-based radio telescopes spanning baseline distances 
from hundreds of kilometres to the Earth’s diameter (see Methods 
for a description of the array). The highly eccentric orbit of the space 
radio telescope, with an apogee of -350,000 km, provided us with 


ground-space fringe detections of the source up to a projected baseline 
distance of eight Earth diameters, probing a wide range of spatial fre- 
quencies perpendicular to the jet propagation direction (Extended Data 
Fig. 1). At the longest projected baselines to RadioAstron, we achieved a 
resolving power of 27 pas, similar to that obtained by the Event Horizon 
Telescope (EHT) at 1.3 mm (-20 pas)’. The large number of detections 
reported within the RadioAstron active galactic nuclei (AGN) survey 
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Fig. 1| The filamentary structure of the jet in3C 279 revealed by RadioAstron. 
a, Total intensity (left) and linearly polarized (right) RadioAstron image at 1.3 cm 
obtained on 10 March 2014. While both images in a show brightness temperature 
(colour scale), the image on the right also shows the recovered electric vector 
position angles overlaid as ticks. Their length and colour are proportional to the 
level of linearly polarized intensity and fractional polarization, respectively. 

b, The 1:1scale 1.3 mm EHT image obtained in April 2017. Contours correspond to 
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our RadioAstron image, and are shown to compare the different scales probed. 
These start at 90% of the peak brightness and decrease by successive factors 

of 3/2 until they reach 5%. Both images were aligned with respect to the pixel 

with maximum brightness. c, The 7 mm VLBA-BU-BLAZAR programme image 
obtained on 25 February 2014. White ellipses at the bottom-left corners of b and c 
indicate the 20 x 20 pas? and 150 x 360 pas? convolving beams, respectively. 

The colour bars refer only to information displayed in a. 


programme’ made 3C 279 an ideal target for detailed imaging. Figure 1 
presents our RadioAstron space VLBI polarimetric image of the blazar 
3C 279. A representative image reconstruction obtained using new 
regularized maximum likelihood methods*” is shown along with the 
closest-in-time 7 mm VLBA-BU-BLAZAR programme image obtained on 
25 February 2014, and the 1.3 mm EHT image obtained in April 2017. We 
show afield of view of around 1 x 1 mas? with an image total flux density 
of 27.16 Jy, and note that all extended emission outside this region is 
resolved out by RadioAstron. The robustness of our image is demon- 
strated in Extended Data Fig. 2, where we show howit fits the data used 
for both total intensity and linearly polarized image reconstruction. 
We acknowledge, however, that VLBI imaging is an ill-posed problem, 
and any image reconstruction that fits the data is not unique (for exam- 
ple, see the comprehensive image analysis carried out in ref. 6). The 
image in Fig. 1is complemented by the 48 images presented in Extended 
Data Fig. 3. In Extended Data Fig. 4 we also demonstrate our ability 
to reconstruct thin filamentary structures with RadioAstron using 
synthetic data (see the Methods for a detailed description). 
Incontrast to the contemporaneous 7 mmand classical millimetre- 
and centimetre-wave VLBI jet images”*, in which the observed syn- 
chrotron emission seems to be contained in a funnel with a uniform 
cross-section, we showin great detail the internal structure of a blazar 


jet and find strong evidence for the filamentary nature of the emitting 
regions within it. We identify the jet core as the upstream bright com- 
ponent, and the so-called core region encompasses approximately the 
inner 200 pas, roughly the extent of the features probed at 1.3 mm. The 
base is slightly elongated and tilted inthe southeast-northwest direc- 
tion, as reported in ref. 2. However, contrary to the EHT sparse sampling 
of the Fourier plane”’, the ground array supporting our space VLBI 
observations provided a substantially larger filling fraction, which ena- 
bled us to reconstruct images with a dynamic range that is two orders 
of magnitude larger. Thus, while we can recover up to three different 
filaments emanating perpendicularly from the jet base, the EHT could 
only recover one and is ‘blind’ to the extended, filamentary structure, 
primarily due to the lack of short baselines. In fact, the peak flux density 
ratio between the core and extended jet measured at 7 mm changes 
from -1.5 in February 2014 to 210 in April 2017. If aligned with respect 
to the brightness peak, both images match remarkably well, and the jet 
feature observed at 1.3 mm is coincident in position and extension with 
our central filament, ignoring the small (<35 pas) core shift between 
the two frequencies”. Within our uncertainty, we do not measure a 
substantial change in the core position angle with respect to the EHT 
image, taken 3 years later. The single-epoch results presented here do 
not allow us to discern whether this elongated structure corresponds 
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Fig. 2 | Analysis of the recovered filamentary structure. a, Left: de-projected 
filament coordinates fitted using three Gaussian curves (darker shading; 
Methods). Lighter shading indicates regions where the coordinates have been 
interpolated. The marker size scales with the flux density. Note that the x- and 
y-axis scales are different. Right: position of the fitted filaments overlaid on the 
reconstructed image. Colour scale as in Fig. 1. The black dotted curve denotes the 
main jet axis. b, Doppler factor computed for plasma propagating along a double 
thread, as is expected for elliptical modes, with / = 13 and 0 =1.9°. 


to the accretion disk or to another extended jet component. 
Nonetheless, based on the small viewing angle inferred” (0 = 1.9°) and 
the multi-epoch kinematic analysis of the model-fitted jet components, 
ref. 2 raised the possibility of this structure corresponding toa sharply 
bent part of the inner jet. 

Moving beyond the core region, in Fig. 2a we show the de-projected 
and on-sky coordinates of the two main filaments (hereafter g and b), 
and a possible third filament (r), obtained from the fitting of three 
Gaussian curves to transverse cuts to the main jet axis (Methods). Fur- 
ther downstream, filament g is continuously recovered and contains 
most of the eastern extended structure flux density. Initially propagat- 
ing in the southern direction, it displays a sharp bend of ~45° to the 
west, close to the core region boundary. Although not continuously, 
we are also able to reconstruct filament b beyond the inner 200 pas 
in what seems to be a helical-like morphology. These two filaments 
converge at ~500 pas down the jet, where filament b crosses over g. 
Further downstream, they bend and converge again at -950 pas, where 
the brightness of the weaker filament is largely enhanced as it bends, 
now dominating the reconstructed emission in the southernmost jet 
region. Some diffuse emission is also systematically recovered parallel 
to filament g after the first crossing, which might indicate the presence 
of athird filament (r). 


As detailed in the Methods, there is a physically consistent domain 
in the space of parameters (F, a,/a.,; where a;is the sound speed of the 
jet flow and a,, is the sound speed of the ambient medium), where 
the Lorentz factor F is within the previously determined range of jet 
flow Lorentz factor for 3C 279 (ref. 12) (r e [10, 40]), which allowed us 
to interpret the observed approximate spatial periodicity, Am, as the 
wavelength of the elliptical surface mode of a kinetically dominated, 
cold jet. More intriguing is the fact that the filaments associated with 
the elliptical mode are brighter in particular locations separated by 
half a wavelength (~400 pas for filament g, ~850 pas for filament b), 
just before the crossings of the filaments. The properties of the flow 
(for example, pressure, density, flow velocity) are modified locally 
by the elliptical wave, with the magnitude of the changes depending 
on the position and time as modulated by the wave phase. Such small 
changes in the properties of the flow could explain the differences in 
brightness between regions inside the jet and, in particular, along the 
filaments. Here the perturbation in the three-velocity vector and the 
subsequent changes in the local Doppler boosting play a major role. 
Figure 2b shows how the Doppler factor, ô, evolves along two threads 
originated by an elliptical mode ina plasma characterized by [= 13 and 
0 =1.9°. With a ratio 6,,,,/Omin = 1.7, the brightness of certain regions 
can increase by a factor of ~5, in agreement with the brightness excess 
observed at -400 pas in filament g, and at ~850 pas in filament b. This 
enhanced emission would then be the result of the Doppler boosted 
emission along the line of sight. This interpretation is supported by the 
fact that both enhanced emission regions are found at approximately 
the same phase of the corresponding helical filament, with the local 
flow velocity pointing along the same direction (the line of sight). 

Continuing with this interpretation, it is important to note that 
the enhanced emission regions will not be steady but will propagate 
downstream at a (pattern) speed equal to the wave’s phase velocity. 
The fact that these brighter regions in filaments g and b match the 
jet features observed at 7 mm (components C32 and C31in refs. 11,13, 
respectively) suggests the appealing possibility that these features 
correspond to the propagation of an elliptical perturbation mode, 
and not to the propagation of a shock, as proposed by the standard 
shock-in-jet model”. In particular, the estimated apparent speed of 
~7.23 c (component C31, ref. 11), corresponding to a propagation speed 
of ~0.996 c (T ~ 11) and close to the lower limit of the 3C 279 jet Lorentz 
factor estimates, supports this possibility as these waves propagate 
downstream along the jet with velocities smaller than or equal to that 
of the jet bulk flow (see refs. 15,16 and the Methods). We illustrate our 
proposed model for the jet variability in 3C 279 in the schematic in 
Fig. 3. Likewise, filamentary structures triggered by plasmainstabilities 
could explain the variability observed at radio wavelengths in other 
blazar sources and, insome cases, they could coexist with components 
originating from shock waves. An important implication of our study 
is that the main emission occurs in thin filaments that cover only a 
fraction of the jet cross-section. The timescale of variability of such 
features can therefore be much shorter than the light-travel time across 
the entire jet width, a fact that can help to explain the extremely rapid 
variability reported in gamma-rays”. We note that while our findings 
and proposed model come from the first spatially resolved image of a 
blazar internal structure, similar models based on differential Doppler 
boosting effects have also been suggested from indirect measurements 
(for example, see ref. 18). Inthe case of 3C 279, magnetic reconnection 
has recently been proposed” as the mechanism behind gamma-ray 
flares, as opposed to shock waves or geometric effects. Nonetheless, 
the variability we model, as it originates from the formation of filamen- 
tary structures triggered in a kinetically dominated jet flow, concerns 
only that observed at radio wavelengths beyond the inner core region. 

Theanalysis of the linearly polarized emission captured by RadioAs- 
tron andthe supporting ground array reveals clear signatures of a toroi- 
dal magnetic field threaded to the relativistic jet. The source is mildly 
polarized, with an integrated degree of linear polarization of -10%. 


Nature Astronomy 


Article 


https://doi.org/10.1038/s41550-023-02105-7 


Black 
hole 
Parsec-scale jet 


Jet funnel 


Helical filaments 


Accretion 
disk 
Region of enhanced 
emissivity (component C32) 


C31 


Fig. 3 | The proposed model for the internal jet structure in3C 279. The 
development of plasma instabilities within the jet flow leads to filamentary 
structures with heterogeneous emissivity. Local changes in the plasma 
properties, which evolve and travel down the jet with the wave velocity, can 
explain these commonly observed moving components as a result of differential 
Doppler boosting. 


Theelectric vector position angle indicates a magnetic field predomi- 
nantly perpendicular to the flow propagation direction—thatis, consist- 
ent witha helical magnetic field dominated by its toroidal component. 
However, the presence of a helical magnetic field usually implies that 
Faraday rotation measure gradients are found across the jet width”? ”. 
Future joint analysis of RadioAstron 22 GHz and VLBA 43 GHz polari- 
metric images will further test this. Relativistic magnetohydrodynamic 
simulations of jets at parsec scales have shown that, inthe presence of a 
helical magnetic field, the observed synchrotron emission is unevenly 
distributed across the jet width” *°. Based on the strong asymmetry in 
the reconstructed emission between the eastern and western sides of 
the jet axis, and following the analysis described in the Methods, we 
can infer a jet bulk flow Lorentz factor of F ~ 13, which is in excellent 
agreement with the estimates provided by analysing the kinematics of 
the parsec-scale jet. Moreover, this allows us to further infer a possible 
helical magnetic field, with an intrinsic pitch angle of ~45°, rotating 
clockwise as seen in the direction of flow motion. 

The findings presented in this Article, supported by previous VSOP 
(for example, refs. 27,28) and RadioAstron (for example, refs. 29,30) 
space VLBI observations, strongly suggest that blazar jets have acom- 
plex and rich internal structure beyond the funnel-like morphologies 
reported by ground-based VLBI studies at lower angular resolutions. 
Future space VLBI missions and enhanced millimetre-wave global 
arrays, which will enable high-dynamic-range observations capable of 
spatially resolving the jet width, should prove decisive in determining 
the true nature of jets powered by supermassive black holes. 


Methods 

Observations 

Observations of 3C 279 were conducted at 22.2 GHz (1.3 cm) on 2014 
March 10-11, spanning a total of 11 h 44 min from 14:15 to 01:59 UT. 
During the observing session, RadioAstron recorded evenly spaced 
(every 80-90 min) blocks of data of 30 min and one final block of 
~2 h, corresponding to its orbit perigee. This allowed the spacecraft to 
cool down its high-gain antenna drive between observing segments. 
Together with RadioAstron, a ground array of 23 antennas observed 
the target: ATCA, Ceduna, Hobart, the Korean VLBI Network antennas 
Tanman, Ulsan and Yonsei, Mopra, Parkes, Sheshan, Badary, Urumqi, 
Hartebeesthoek, Kalyazin, Metsadhovi, Noto, Torun, Medicina, Onsala, 
Yebes, Jodrell Bank, Effelsberg, Svetloe and Zelenchukskaya. 

Left- and right-circularly polarized (LCP and RCP, respectively) 
signals were recorded simultaneously at each station, with a total 
bandwidth of 32 MHz per polarization. The collected data were then 
processed at the Max-Planck-Institut fiir Radiostronomie using the 
upgraded version of the DiFX correlator”. Fringes between Radio- 
Astron and ground stations were searched using the largest dishes, 


separately for each scan. This provided a first-order clock correction 
that was later refined with baseline stacking in AIPS”. When no signal 
was found, we adopted a best-guess clock value extrapolated from 
scans giving fringes, with the aim of performing a further global fringe 
search at a later stage with AIPS. 


Data reduction 

For the initial data reduction, we used ParselTongue®, a Python inter- 
face for AIPS. Inthe first stage, we performed ana priori calibration of 
the correlated visibility amplitudes using the system temperatures and 
gain curves registered at each station. Some of the antennas participat- 
ing in the observations failed to deliver system temperature informa- 
tion, which we compensated by using nominal values modulated by 
the antenna’s elevation at each scan. As we chose the average system 
temperature as the station’s default value, visibility amplitudes were 
not properly scaled. We overcame this issue by determining, every 
2 min, the gain corrections needed for each intermediate frequency and 
polarization froma preliminary image in which only closure quantities 
(closure phases and log-closure amplitudes) were involved, using the 
SMILI software library****. We then applied to each antenna the mean 
gain value obtained, allowing further residual corrections to be made 
during the final imaging and self-calibration. The image total flux 
density was fixed to that measured by the intra-Korean VLBI Network 
baselines (27.65 Jy), whose a priori calibration was excellent”®. Finally, 
we corrected the phase rotation introduced by the receiving systems 
as the source’s parallactic angle changes. 

We then solved for residual single- and multi-band delays, phases 
and phase rates by incrementally fringe-fitting the data. In the first 
iteration we excluded RadioAstron and performed a global fringe 
search on the ground array witha solution interval of 60 s using Mopra 
and Effelsberg as reference antennas for the first and second part of 
the experiment, respectively. Once fully calibrated, the ground array 
was coherently combined (trough baseline stacking) to increase the 
signal-to-noise ratio of possible fringe detections to RadioAstron. 
To account for the acceleration of the spacecraft near its perigee and 
the low sensitivity of the longest projected baseline lengths to it, we 
adopted different solution intervals (from 10 s to 240 s) and data 
total bandwidth (by combining intermediate frequencies). With a 
signal-to-noise ratio cutoff of 5, reliable ground-space fringes were 
detected at distances of up to ~8 Earth diameters, corresponding tothe 
first observing block of RadioAstron (from 14:15 to 14:45 UT), providing 
amaximum angular resolution of 27 pas in the direction transverse to 
the jet axis. Finally, we solved for the antennas’ bandpasses, the delay 
differences between polarizations using the task RLDLY, and exported 
the frequency-averaged data along each intermediate frequency. 
The fringe-fitted visibility coverage in the Fourier plane is shown in 
Extended Data Fig. 1. 


Imaging 

Imaging of the data was carried out using novel regularized maximum 
likelihood (RML) methods”, implemented in the eht-imaging software 
library’. While the CLEAN algorithm” has been widely used in the 
past for VLBI image reconstruction, RML methods are not extensively 
used, especially at centimetre wavelengths and for space VLBI experi- 
ments. Generally speaking, RML methods try to solve for the image 
/, constrained by the measurements V, that minimizes an objective 
function/ defined as follows: 


D= 5 


data terms 


anxo V- >) BrSROD, a) 


reg. terms 


where D and Rare a set of selected data products and regularization 
terms, respectively, and a, and 6, are hyperparameters that weight the 
contribution of the image fitting to the data y”,, and the image-domain 
regularization S,, to the minimization of the previous equation. 
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Contrary to traditional CLEAN algorithm, full closure data products 
(closure phases and log closure amplitudes) can be employed during 
image reconstruction in addition to complex visibilities, further con- 
straining the proposed image. Given the large number of telescopes 
participating in the experiment, closure quantities have proved quite 
useful as atmospheric phase corruption and gain uncertainties are 
mitigated. Multiple regularization over the proposed image can be 
imposed too, such as smoothness between adjacent pixels or similarity 
toaprior image. 

Before imaging, we performed an initial phase-only self-calibration 
to a point source model with a solution interval of 5 s and coherently 
averaged the data in 120 s intervals using the DIFMAP package”. We 
compared these results with those obtained with the AIPS task CALIB, 
for which a signal-to-noise ratio cutoff of 5 was set, to ensure that no 
artificial signal was introduced into the data. In the following para- 
graphs we describe the imaging procedure. 

Asafirsta step, we flagged all baselines to RadioAstron and imaged 
the data collected only by ground radio telescopes. The pre-processed 
data noise budget was inflated by a small amount (1.5%), to account 
for non-closing errors, and the image was initialized with an ellipti- 
cal Gaussian oriented at roughly the same angle as the 7 mm image 
and enclosed in a 1.5 x 1.5 mas? field of view gridded by 200 x 200 
pixels. As mentioned above, because of the poor a priori amplitude 
calibration due to missing antennas’ system temperatures, we opted 
for a first round of imaging in which only closure quantities (closure 
phases and log closure amplitudes) were used to constrain the image 
likelihood. This likelihood took the form of the mean squared standard- 
ized residual (similar to a reduced x°) as defined in ref. 6. Each imaging 
iteration took the image reconstructed in the previous step as its initial 
guess blurred to the nominal resolution of the ground array (that is, 
223 pas), which prevented the algorithm getting caught in local minima 
during the optimization of equation (1). We then self-calibrated the 
data to the closure-only image obtained and incorporated full complex 
visibilities in the imaging process, which was finalized by repeating 
the imaging and self-calibration cycle two more times. In addition to 
the data products mentioned, we imposed several regularizations 
on the proposed images. These included maximum entropy (mem), 
which favoured similarity to a previous image; total variation (tv) and 
total squared variation (tv2), which favoured smoothness between 
adjacent pixels; 2, - norm, which favoured sparsity in the image; and 
total flux regularization, which encouraged a certain total flux density 
in the image. Finally, we restored all baselines to RadioAstron and 
repeated this procedure, substituting the Gaussian initialization with 
the blurred, ground-only image previously reconstructed and using 
the full array nominal resolution (27 pas) between imaging iterations 
to blur intermediate reconstructions. 

In contrast to full Bayesian methods, RML techniques do not esti- 
mate the posterior distribution of the underlying image, but instead 
compute the maximum a posteriori solution; that is, the single image 
that best minimized equation (1). The hyperparameters chosen will 
inevitably have an impact on the reconstructed image features, thus 
we conducted a scripted parameter survey to ensure the robustness of 
the subtle structures seen in Fig. 1and to impartially determine which 
parameters performed better in the image reconstruction. From the 
many images obtained, we show in Extended Data Fig. 3 the complete 
collection of images that could potentially describe the observed 
source structure. These fitted the data equally well and preserved 
the total flux measured by the Korean VLBI Network to a certain level. 
The regularizers and hyperparameters used to obtain these images 
are listed in each panel of Extended Data Fig. 3. Apart from these, we 
gave the same weight to complex visibilities and closure quantity data 
terms. Although we observed some differences in the weaker emission, 
the main filamentary structure was present in all the images. Figure 
lcorresponds to no. 21, which had the overall minimum reduced y’. 


Synthetic data tests 

Given that one of the helical filaments reconstructed in our images 
(filament b) crosses over the other with a position angle roughly ori- 
ented in the direction where we lack space-ground baseline cover- 
age, we tested the ability of the image reconstruction algorithm to 
successfully recover similarly oriented structures. With that aim, we 
generated two synthetic datasets using eht-imaging (for example, see 
refs. 6,40). These simulated space VLBI observations of the source 
models presented in Extended Data Fig. 4 with the same baseline cover- 
age as that in Extended Data Fig. 1. The datasets were corrupted with 
thermal noise and the visibility phases were scrambled. These simple 
geometric models were specially designed to mimic the crossing of two 
helical filaments and to some extent the brightness distribution of the 
3C 279 images reconstructed. While the first model was oriented as the 
parsec-scale jet in3C 279, the second model was oriented perpendicu- 
larly to the fitted beam position angle (~10°). 

In Extended Data Fig. 4 we show the images reconstructed 
following the procedure described in the previous section. For 
comparison, we also present the images obtained after flagging all 
RadioAstron baselines; that is, with only the ground array. When the 
full array was considered, we were able to fully recover the crossing 
of the filaments and ground-truth structure of the model oriented 
in the same way as the jet in 3C 279, although we acknowledge the 
poorer performance in reconstructing the internal structure of 
the model oriented perpendicularly to the position angle of the 
fitted beam, as expected. On the contrary, when space-ground 
baselines are removed, the images reconstructed display a jet cor- 
rectly oriented but with no traces of helical filaments. Thus, we are 
confident that the recovered filamentary structure, only achievable 
at this observing frequency thanks to RadioAstron, is derived from 
the intrinsic source structure and not from the lack of north-south 
space-ground baselines. 


Polarimetric imaging 

The polarization results presented in Fig. 1 were also obtained using 
the eht-imaging library. A more complete description of the method 
can be found in refs. 4,41, here we briefly outline the procedure fol- 
lowed. For polarimetric imaging, eht-imaging minimizes equation (1), 
substituting complex visibilities and closure quantities data terms by 
polarimetric visibilities ? = 6 + i ŭ and the visibility domain polari- 
metric ratio m = ?/J, constructed using the visibility-domain Stokes 
parameters /, Q, and U. Note that total intensity and linearly polarized 
intensity images were reconstructed independently. Image regulariza- 
tion then included total variation, which, as for total intensity imaging, 
encouraged smoothness between adjacent pixels, and the Holdaway- 
Wardle regularizer®, which prefers pixels with polarization fraction 
values below the theoretical maximum of 0.75. The pipeline then alter- 
nated between minimizing the polarimetric objective function and 
solving for the complex instrumental polarization, the so-called D 
terms. The instrumental polarization calibration was performed by 
maximizing the consistency between the self-calibrated data and 
sampled data from corrupted image reconstructions. After D-term 
solutions were found for each antenna, the reconstructed polarimetric 
image was blurred, as was done for Stokes 7 imaging, and the imaging- 
calibration cycle was repeated until the solutions converged. 

In addition to instrumental polarization, VLBI polarimetric analy- 
ses rely onthe calibration of the absolute polarization angle. To account 
for this, we compared our polarization results with the closest-in-time 
7 mm results. The recovered polarization angle patterns matched 
remarkably well when our image was convolved with the 7 mm beam. 
Ignoring Faraday rotation of the polarization angle between the two 
frequencies, based on the small rotation measure values reported in 
refs. 43,44, we estimated an overall median difference of ~8° that we 
applied to the results presented in Fig. 1. 
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Filament fitting 

The relative right ascension and declination of the filaments were 
obtained by fitting three Gaussian curves to transverse profiles of 
the brightness distribution. We first computed the main jet axis, 
commonly referred to as the ridge line, from a convolved version of 
the reconstructed image. Using a sufficiently large Gaussian kernel, 
we blurred our image until the emission blended into a unique stream 
and the filaments were no longer distinguishable, similarly to the 
7 mm VLBA-BU-BLAZAR image. We then projected this image into 
polar coordinates, centred at the jet origin, and sliced it horizontally, 
storing the position of the flux density peak for each cut. These posi- 
tions were then transformed back to Cartesian coordinates, provid- 
ing the main jet axis. To each pair of consecutive points conforming 
the axis, we computed the local perpendicular line and retrieved the 
flux density of the pixels contained in the cut. With this procedure, 
we assembled a set of transverse brightness profiles to which we 
fitted the sum of three Gaussian curves using the python package 
Imfit*’. The number of Gaussian components used was motivated by 
the number of filaments observed emanating from the core region, 
although we note that two Gaussian components were enough to 
fit the two main threads. Finally, we selected the coordinates as the 
position of the peak(s) found in the curve best fitting each cut. In 
Fig. 2, coordinates were de-projected assuming a source redshift 
z= 0.536 (ref. 46), a viewing angle 0 =1.9° (ref. 11) and a cosmology 
H, = 67.7 km s” Mpc”, Q,, = 0.307 and Q, = 0.693 (ref. 47), where Hois 
the Hubble constant, and Q,, and Q, are the matter and dark energy 
density parameters, respectively. 


Instability analysis 

Based on the aforementioned Gaussian fitting to the observed fila- 
ments, we estimated an approximate spatial periodicity A,, of 950 pas 
(projected on the plane of the sky) or 175 pc (de-projected), which 
corresponded to -2.3 x 10° gravitational radii assuming a black hole 
mass of Mpu ~ 8 x 10° M, (ref. 48). The possibility of these filaments 
reflecting a fundamental periodicity of the black hole or inner accre- 
tion disk directly associated with their rotation should be dismissed, 
as it would imply propagation speeds along the filaments larger than 
the speed of light by orders of magnitude. At the same time, explana- 
tions of such a fundamental periodicity in terms of the precession of 
ajet nozzle, caused by the Lense-Thirring effect*’ or a supermassive 
black hole binary system, invoked to explain a sharp bend inthe nuclear 
region, have recently been refuted’. On the other hand, anchoring the 
filaments to the outer accretion disk to allow a subluminal propagation 
of the helical pattern would imply an exceedingly large (Keplerian) 
disk radius—that is, larger than ~1 light yr, which is about two orders 
of magnitude larger than the expected disk sizes®. 

According to ref. 2, the jet no longer accelerates beyond ~100 pas 
from the core, suggesting a kinetically dominated flow in which the 
observed filaments show a magnetic field structure dominated by the 
toroidal component. Taking this into account, we conclude that these 
bright filaments reveal compressed regions with enhanced gas and 
magnetic pressure—favouring an increased synchrotron emissivity and 
ordering of the magnetic field. Thus, these might be associated with the 
triggering and development of flow instabilities. Current-driven kink or 
Kelvin-Helmholtz instabilities are the most plausible mechanisms for 
developing such helical structures", Rayleigh-Taylor instabilities 
have also been discussed in the context of jet expansion and recollima- 
tion asa possible trigger of small-scale distortions of the jet surface and 
turbulent mixing (see ref. 52 for a review). However, Rayleigh-Taylor 
instabilities would not produce filaments like the ones we observe, 
so we neglect this option. Current-driven instabilities dominate in 
Poynting-flux regimes with strong helical magnetic fields; thatis, inthe 
jet’s acceleration and collimation region. Kelvin-Helmholtz instabili- 
ties, however, have the largest growth in kinetically dominated flows, 
and are thus favoured in our case. The extension of the filaments greatly 


exceeds the jet radius, which is expected for Kelvin-Helmholtz surface 
modes. While two filaments could be generated by an elliptical mode, 
the possible third filament observed might indicate the presence of an 
additional helical mode interfering with the elliptical. 

Assuming that the jet is kinetically dominated and cold, 
as expected for powerful jets that are already expanded and acceler- 
ated, the fastest growing frequency of a mode is given by 
O% )R/dex = (n + 2m + 1/2)m/2 (refs. 15,53; quantities with a superscript 
asterisk refer to resonant, or fastest growing, modes), where R is the 
radius and nand mare the types of mode (n = 1, 2 for helical and ellipti- 
cal modes, respectively; and m = 0 for a surface mode). Taking 
@ < 2m1c/A,,, we found ae = 10° c for both the helical and elliptical 
modes. At this maximum growth frequency, and for a highly supersonic 
jet (thatis, with a jet Mach number M, > 1), the wavelength of the mode 
and wave velocity are given, respectively, by: 


ee r 
m nti “a/dx +0 


(2) 


r 
vs x —— l, (3) 
“O Gldex +T 


where u is the jet flow velocity (which approximates the light speed c 
given the large Lorentz factors inferred), lis the jet flow Lorentz factor 
and M.x is the Mach number of the jet with respect to the ambient sound 
speed (Mex = U/dex = C/Aex = 100). 

In our interpretation, described in the main text, the wave velocity 
vi coincides with the (pattern) speed of the jet feature observed at 
7mm (component C31in ref. 11), very close toc, leading to the condition 
(see equation (3)) a/a,, < or, equivalently, M, > M.,/T. With Mex = 100 
and/ € [10, 40], the last condition implies M; > 1, hence validating our 
assumption of a kinetically dominated flow at the observed scales. 


Jet properties derived from the reconstructed polarimetric 
emission 

The synchrotron radiation coefficients are a function of the angle 
between the magnetic field and the line of sight in the fluid frame. 
Thus, for a fixed viewing angle and jet flow velocity, the bulk of the 
emission will be located on either side of the main jet axis, depending 
onthe magnetic field helical pitch angle. This asymmetry is maximized 
when the helical magnetic field pitch angle (in the fluid frame) @’ is 45° 
(ref. 23). Given the strong asymmetry in the reconstructed emission 
between the eastern and western sides of the jet axis, we could assume 
that the viewing angle in the fluid’s frame approximates @’ (that is, 
0’ ~ d’). Hence, given the estimated viewing angle in the observer’s 
frame (0 ~1.9°) (ref. 11) and the light aberration transformations” 


sin 0 
Td —Bcos6)’ 1- gcos’ 


sing’ 


(4) 


where £ = y1 — 1/T2,we can infer a jet bulk flow Lorentz factor of [= 13, 
whichis in excellent agreement with the estimates provided by analys- 
ing the kinematics of the parsec-scale jet", and satisfies the upper limit 
previously established by our Kelvin-Helmholtz instability analysis. 
Moreover, this allowed us to estimate the viewing angle 6, at which 
the emission asymmetry will reverse from one side to the other as 
cos 6, = (1— yr)” (ref. 23), which results in 0, ~ 4.4° for P= 13. Given 
that 0 < 6, and the bulk of the reconstructed emission is located to the 
east of the jet axis, we infer a helical magnetic field rotating clockwise 
as seen in the direction of flow motion. The Lorentz transformation of 
the magnetic field from the fluid’s to the observer’s frame boosts the 
toroidal component by T (ref. 55), and therefore the helix pitch angle 
transforms as tan @ = T tan ġ'. This means that @ ~ 86° in the observer 
frame, whichis in agreement with the predominantly toroidal magnetic 
field observed. 
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Data availability 
The pre-processed dataset used for imaging is available via Github at 
https://github.com/aefez/radioastron-3c279-2014. 


Code availability 

The software packages used to calibrate, image and analyse the data 
are available at the following websites: AIPS, http://www.aips.nrao. 
edu/index.shtml; ParselTongue (https://www.jive.eu/jivewiki/doku. 
php?id=parseltongue:parseltongue), DIFMAP (https://science.nrao. 
edu/facilities/vlba/docs/manuals/oss2013a/post-processing-software/ 
difmap), SMILI (https://github.com/astrosmili/smili), eht-imaging 
(https://github.com/achael/eht-imaging) and Imfit (https://Imfit. 
github.io/Imfit-py/). 
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Extended Data Fig. 1| Baseline coverage for our RadioAstron observations of 3C 279 in March 2014. Rainbow-coloured and grey points indicate individual ground- 
ground baselines and space-ground baselines, respectively. Dashed circles indicate the baseline length in Earth’s diameter units (Dẹ) and the corresponding angular 


resolution. 
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Extended Data Fig. 2 | Fitting of the polarimetric RadioAstron image to a function of time. Error bars indicate + lo uncertainty from thermal noise plus 
selection of data products. Two minute time-averaged data (black points) 1.5 % non-closing error uncertainties added in quadrature. All these examples 
and image model (red points) self-calibrated visibility amplitudes and phases, include RadioAstron measurements. 
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Extended Data Fig. 3 | Top 48 image reconstructions from the parameter survey conducted. Each image includes the closure phase (cp) and log closure amplitude 
(Ica) reduced x’, the image regularizers used and their weight, and the total flux reconstructed. 
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Extended Data Fig. 4 | Synthetic data tests. In the top row we present the geometric models used to generate the synthetic data. In the middle and bottom rows we 
show, respectively, the images reconstructed from each data set when RadioAstron is included in the array and when only ground stations participate. 
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